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ABSTRACT 


A ship experiences different sea conditions during sailing. Particularly in 
adverse conditions, the ship generator will be put under extreme loading stress due 
to power fluctuations. This could potentially damage or shut down the generator. 
Having a Hybrid Energy Storage System (HESS) reduces the loading stress on the 
generator, as it can take over as an alternative power source. 

The goal of the thesis is to design a control strategy using Model Predictive 
Control (MPC) to effectively control the HESS output behaviors in order to 
compensate for power fluctuations due to ocean waves. For the purpose of this thesis, 
the HESS is composed of a combination of batteries and supercapacitors. 

Computer simulation using MATLAB verified that the MPC is able to 
effectively control the HESS in compensating for power fluctuations and at the same 
time, minimize generator overload. Different weightages on the battery were also 
tested to observe the battery and supercapacitor output behaviors. The simulation also 
provided insight into the impact of the weightages on the states and inputs, and it was 
observed that there is saturation on the battery and supercapacitor currents when their 
weightage difference is greater than 20. In addition, the model could also be used 
for preliminary study for HESS sizing based on the project power demands. 
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I. 


INTRODUCTION 


According to the United Nations Conference on Trade and Development 
(UNCTAD), the maritime industries account for 70% of world trade [l]-[3], and based on 
a 2012 study, the shipping industry produced about 796 million tons of carbon dioxide 
(C02) [4], [5]. Based on the projected economic growth and energy development of these 
industries, the International Maritime Organization (IMO) forecasts that C02 emissions 
could increase from 50% to 250% by 2050 [4]. The IMO states that to further reduce these 
industries environmental footprint, the IMO has proposed reducing their emissions from 
3.5% to 0.5% by January 2020. It has also set targets for the shipping industry to cut down 
Green House Gas (GHG) emissions by at least 50% by 2050. 

Technologies such as hybrid energy storage systems (HESS), variable frequency 
drives, and demand-side management, which are widely used in land-based microgrids and 
hybrid vehicles, are starting to emerge in the marine industries. While a high percentage of 
sea-going vessels are still diesel-powered, there is an increasing trend of ship owners and 
operators turning to greener technologies that can help them in reducing fuel bills and 
environmental impact. In particular, electrification in the marine industry has become an 
attractive technology, and the U.S Navy [5] has proposed the concept of an All Electric 
Ship (AES). Several military agencies have also embarked on AES research and designs 
[6]-[8]. Applying HESS in ships can generate potential benefits such as increasing 
flexibility in ship design, addressing power fluctuations, improving generator reliability, 
increasing fuel efficiency, and decreasing the C02 footprint. 

A hybrid energy storage system is defined as a combination of two or more types 
of power sources/energy storage devices (ESD) that can provide energy [9]. Generally, 
they are classified into three categories: 

1. Traditional internal combustion engines (ICE) such as diesel generators, 
gas, or steam turbines; 

2. Electrochemical energy storage such as batteries and power sources such 
as fuel cells; 
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3. 


Other energy storage such as flywheels and supercapacitors. 


ESDs can offer many benefits to the marine vessels. As compared to diesel engines, 
ESDs can quickly adapt to sudden changes in load, hence improving the stability of the 
overall powertrain system. Also, ESDs can reduce engine maintenance, optimize fuel 
consumption and minimize the power generation machinery footprint on the ship [10]. A 
previous thesis by Raison [11], also examined the implementation of supercapacitors (SC) 
within the energy storage system of an electric propulsion system (EPS) on a hybrid- 
powered vessel. 

ESDs are also effective in load peak shaving, load leveling, power smoothing, 
handling frequency and voltage fluctuations, and ensuring power quality. Development in 
ESDs also allows users to design a more energy-efficient ship as well as offer flexibility in 
power supply configuration. Unlike electric vehicles, ships can have several independently 
controlled energy-generating sources connected to a standard switchboard. 

Several power management control strategies are commonly applied in the 
industry, including static optimization methods [12], [13], dynamic optimization methods 
[14]-[16], and predictive powertrain controls for minimizing fuel consumption [14]-[15]. 
With proper control strategies, the ESDs can complement each other and perform 
effectively. In [10], a detailed review of the types of smart ships is provided and a 
description of the key benefits and challenges of the different control strategies between 
the propulsion types. 

In hybrid energy systems, the control strategy used is also a key component. A 
suitable control strategy that can forecast the energy demand and compute in real-time the 
optimal energy distribution among each ESD makes the microgrid more energy efficient. 

In recent years, Model Predictive Control (MPC) has been gaining popularity in 
HESS energy management applications [17]. MPC is a model-based predictive control 
method that uses the first reference values and projects them over a prediction horizon 
while minimizing a cost function and satisfying a set of hard constraints on plant behavior. 
This recurring computation is performed based on defined intervals. MPC has been applied 
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in a broad range of applications, such as engine control, energy management, etc. [18]- 
[23], 

This research studies an MPC method that minimizes the change in steady-state 
velocity when navigating through wave as well as minimizing the generator current by 
optimizing the HESS ESD output based on weightage and constraints. The MPC controls 
the optimal energy distribution between the available energy sources over the prediction 
horizon. The proposed algorithm has been simulated using a simplified dynamic model of 
the Norwegian car ferry, MF Oppedal [24]. The simulation results show the different 
tradeoffs between the energy storage devices. 

This thesis is organized as follows: Chapter II introduces the types of ESDs and 
provides an overview of the types of control strategies. Chapter III presents the ship 
parameters calculations as well as the MPC state-space model. Chapter IV discusses the 
scenarios and simulation results. Lastly, Chapter V concludes the research and discusses 
areas for future work. 


3 



THIS PAGE INTENTIONALLY LEFT BLANK 


4 



II. BACKGROUND 


A. ENERGY STORAGE DEVICES 

The operating characteristics of various energy storage devices differ from one 
another in terms of technical features such as power and energy densities, charge and 
discharge time, operating temperature, lifetime, and maintenance requirements. Table 1 
provides typical technical features of various kinds of ESDs [9]. 

Over the years, rapid developments in power electronics have also improved ESD 
capabilities. With the increasing demand for cleaner energy and lower emissions, many 
researchers and original equipment manufacturers are exploring the idea of combining 
traditional internal combustion generators with ESDs for improved efficiency. 

ESDs can be classified into three main categories, i.e., electrical energy storage, 
chemical energy storage, and mechanical energy storage. Figure 1 shows an overview of 
the energy storage devices in these three categories, while Figure 2 depicts the comparison 
of energy density, power density, and energy exchange related to specific ESD types. It is 
observed that batteries have the lowest power density in terms of kWh/kg but the highest 
$/kW This could indicate that a battery would be an expensive option if a large amount of 
power is needed in a microgrid. The following sections will discuss in greater detail the 
battery and SCs operating characteristics and how they complement each other. 


Table 1. The technical features of ESDs. Source: [3]. 


System 

Power Density 

Energy Density 

$/kW 

Efficiency 

Life Time 

Response Time 


(kW/kg) 

(kWh/kg) 


(%) 

(years) 


Lead Acid 

75-300 

30-50 

300-600 

65-80 

3-15 

ms 

NiCd 

150-300 

50-75 

500-1500 

75-85 

5-20 

ms 

NaS 

150-230 

150-240 

1000-3000 

75-90 

10-15 

ms 

Li-ion 

150-315 

75-200 

1200-4000 

90-97 

5-100 

ms—s 

Fuel Cells 

500+ 

800-10,000 

10,000+ 

20-50 

10-30 

ms-min 

SMES 

500-2000 

0.5-5 

200-300 

90-95 

20+ 

ms 

Flywheel 

400-1500 

10-30 

250-350 

90-95 

15-20 

ms-s 

Ultra-capacitor 

100,000+ 

20+ 

100-300 

85-98 

4-12 

ms 
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Figure 1. Energy storage technologies. Source: [9]. 



Figure 2. Ragone performance chart comparing the energy density, power 
density, and duration of energy exchange for different ESDs. Source: [25]. 


1. Batteries 

Batteries transform chemical energy directly into electrical energy through an 
electrochemical oxidation-reduction reaction [9]. They are commonly used in 
Uninterruptible Power Supplies, automobiles, ships, microgrids, and other applications. In 
recent years, the notable increase in demand from the hybrid-electric automotive industry 
has led to greater demands for batteries capable of producing higher energy density. Such 
batteries include the lithium-ion type. Table 2 lists the technical specifications of the 
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lithium-ion and other types of batteries. It can be noted that a lithium-ion battery 
outperforms the other three batteries types, and therefore, the lithium-ion type is also the 
most commonly used. 


Table 2. Comparison of energy delivery profile technologies. Source: [25]. 



Lead-Acid 

Ni-Cd 

Ni-MH 

Li-Ion 

Cell voltage (V) 

2 

1.2 

1.2 

3.2 (LiFeP0 4 )3.6 
(NMC) 

Specific energy (Wh/kg) 

1-60 

20-55 

1-80 

3-100 

Specific power (W/kg) 

<300 

150-300 

<200 

100-1,000 

Energy density (kWh/m 3 ) 

25-60 

25 

70-100 

80-200 

Power density (MW/m 3 ) 

<0.6 

0.125 

1.5-4 

0.4-2 

Maximum cycles 

200-700 

500-1,000 

600-1,000 

3,000 

Discharge time range 

>1 min 

1 min-8h 

>1 min 

10s-l h 

Cost (S/kWh) 

125 

600 

540 

600 

Cost (S/kW) 

200 

600 

1,000 

1,100 

Efficiency (%) 

75 

75 

81 

99 


Table 3 compares the advantages and disadvantages of the different types of 
batteries. Although the lead-acid battery is cheaper, it has a short cycle life and very low 
efficiency as compared to the lithium-ion battery. In addition, the disposal of lead-acid 
batteries at the end of their life cycle presents challenges. 


Table 3. A comparison between battery technologies. Source: [9]. 


Type of Battery 

Advantages 

Disadvantages 

Lead Acid 

Inexpensive 

Lead is easily recyclable 

low self-discharge (2-5% per month) 

Shorty cycle-life (around 1500 cycles) 

Cycle life is affected by depth of charge 

Low energy density (about 30-50 kWh/kg) 

Nickel Cadmium 

High energy density (50-75 kWh/kg) 

High cycle count (1500-3000 cycles) 

High degradation 

High cost 

Toxicity of cadmium metal 

Sodium Sulphur 

High energy density (150-240 kWh/kg) 

No self-discharge 

No degradation for deep charge 

High efficiency (75-90%) 

Temperature of battery is kept between 

300 °C to 350 °C 

Lithium-ion 

Very high efficiency (90-97%) 

Very low self-discharge (1-3% per month) 
Low maintenance 

Very high cost 

Life cycle reduces by deep discharge 

Need special overcharge protection circuit 
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2 . 


Supercapacitors 


Supercapacitors (SC), also known as ultra-capacitors, are devices that can store 
large amounts of electrical charge [9]: 

An SC uses two mechanisms to save electrical energy: double-layer 
capacitance and pseudocapacitance. Double-layer capacitance is 
electrostatic in origin, while pseudocapacitance is electrochemical, which 
means that supercapacitors combine the workings of regular capacitors with 
the operations of an ordinary battery. [27] 

An SC also has a higher energy density than regular capacitors. Furthermore, SCs have low 
impedance, which enables them to supply power efficiently. They are generally applied 
where high energy is required for a shorter time. SCs are often paired with batteries as a 
complementary storage device to prolong the shelf life of the battery. Figure 3 shows the 
individual structure of an SC. 


-AV- 


Current Collector 




Porous electrode 




Electrolyte 




Separator 


- 1 

► 


Figure 3. Supercapacitor cell. Source: [9]. 


Among the key features of SCs are higher power density, faster discharging and 
charging, and enhanced life cycle, but they often entail higher cost/Wh [9]. Other 
disadvantages of SCs include low energy density, high self-discharge, overcharging, and 
low cell voltage. 

The pairing of SCs and batteries is commonly seen in hybrid vehicles [8-12] as the 
combination can improve the lifetime, performance, and lifecycle of the battery. This 
pairing has been extensively proposed in the research literature [3], [15], [16], [20], 
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[31]—[34]. The battery and SC complement each other well based on their contrasting 
energy and power densities. The SC also possesses a much longer lifecycle than the battery. 
This combination is usually considered for most HESS microgrids mainly due to their 
similarity in working principle, relatively low cost, and their complementary characteristics 
[9], [10]. This complementary pairing is used for matching the energy demands studied in 
this research project, which simulates the power demands of a ship navigating through 
different profiles of sea conditions. 

3. Load Leveling and Peak Shaving 

The creation and storage of energy during low power demand periods and the use 
of that energy during peak demands is defined as “load leveling” [35]. This approach is 
useful in reducing energy fluctuations, as shown in Figure 4. By contrast, moderating and 
reducing the peaks is called “peak shaving.” Figure 5 shows an illustration of peak shaving 
in a marine vessel. ESDs, especially batteries, have been widely adopted to cater to variable 
load in hybrid vehicles, and this application has been extensively researched in the 
literature [31], [36]. 



Figure 4. Example of load leveling application. Source: [9]. 
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Figure 5. Peak shaving application in marine vessel. Source: [9]. 


B. MODEL PREDICTIVE CONTROL 

The development of hybrid powertrains has been evolving rapidly in the last two 
decades [37]. HESS powertrains, consisting of two or more energy sources, are required to 
deal with different energy forms to deliver the required power demand. In this kind of 
system, a primary controller coordinates the energy sources so that the power demand can 
be satisfied efficiently and effectively. Efficient power management can lower fuel 
consumption and emission rates of the hybrid powertrain. In general, power management 
control strategies can be classified into rule-based control strategies and optimization- 
based control strategies [38], as illustrated in Figure 6. 
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Power management strategies 


Rule-based (RB) 


Optimization-based 


Deterministic RB Fuzzy RB Global optimization 




1 

1 

Thermostatic 

on/off 

Power 

follower 

n 

Adaptive 

RB 

Ada; 

5 five 

1 1 

Dynamic Static 

optimization optimization 


Frequency- Point/s- 

,. Basic Predictive 

based tracking 


Real-time optimization 

I 


r 

PMP 


RC 


ADP 


MPC ECMS ES 


Figure 6. Classification of power management strategies. Source: [38]. 


1. Rule-Based Strategy 

As explained in [38], rule-based techniques are ad hoc strategies in which the 
control law is defined by on-off or by fuzzy-logic rules [39]. Rule-based strategies are 
directly structured and easily implementable and are usually devised using human know¬ 
how and trial and error. Rule-based strategies are commonly used by hybrid vehicle 
manufacturers as this approach has low computational requirements [40]. More in-depth 
discussion on the different types of rule-based strategies is highlighted in [38]. It is 
important to note here, however, the main drawback of rule-based strategies is that they 
are unable to effectively react to power demands that are not defined by the rules. 

2. Optimization-Based Strategy 

Optimization-based strategies have been gaining more researcher attention and are 
also widely applied in hybrid powertrains [40]. They are based on optimal control 
techniques that minimize a suitable cost function, usually a quadratic function. In 
particular, techniques based on optimization in a finite time prediction horizon, such as the 
Model Predictive Control (MPC), are well known in the literature. They can be 
implemented in real time to provide prompt handling of power demands. For this thesis, as 
the load demand is dynamic and uncertain, MPC is applied as the control strategy for the 
hybrid powertrain. 

MPC is gaining popularity and has been more widely used over the last two 
decades, both within the research control community and in hybrid microgrid industries 
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[41]. The model used for prediction is the standard discrete-time state space model of the 
form [42]: 

x(t+ 1) = Fx(t) + Gu(t) (1) 

y(t ) = Hx(t ) + Ju(t) (2) 

Hence, at any time t and for the given state vector x(t) and cost function 
g(x(t),x l+Ut ,...,x l+Nlt u Ilt ,...,u t+m ) associated with predicted states x t+klt =x(t + k), k=l..,N, 
and inputs u t+k]t = u(t + k ) , k=0,..,N, can be defined. This cost function, usually quadratic, 
can then be optimized with respect to the control sequence u Hk[t , from which we determine 
the control u(t) = u t]t to be applied at the present time t. 

The MPC principle is illustrated in Figure 7. As it is shown, the control input u(t) 
at time t is based on a finite time horizon prediction based on the current state x(t) and 
optimized with respect to future inputs u t+a . 


Model Predictive Control 


Objectives Mode) Constraints 




Reference 

Optimizer 

B 

Input 


** nani 


1—=J 



A 


Measurements 


Output 



Receding horizon strategy introduces feedback 


Figure 7. An illustration of MPC working principle. Source: [43]. 
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Furthermore, MPC will also consider system limitations simply by establishing 
them as constraints in the optimization problem. Examples of these limitations could be 
the battery state of charge or the maximum charge power available in a battery [44]. 

For this thesis, the load model for the hybrid ship is dynamic and is supported by 
three energy sources. MPC enables us to determine the energy sources contribution 
weightage as well as to set limitations on each energy source. 
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III. PROPOSED ARCHITECTURE 
AND MATHEMATICAL MODELS 


This chapter introduces the dynamic characteristics of the various components and 
the proposed architecture for this research. In particular, it focuses on the dynamics of a 
particular ship, the MF Oppedal, which we will be using to test the effectiveness of the 
proposed algorithm. In addition, we will be introducing dynamic models of relevant 
components, such as batteries and generators, of the Hybrid Energy System. The goal is to 
formulate a set of state space equations and constraints, which will be used to design the 
optimization algorithm. 

A. SHIP MODEL AND SHIP DRAG (RESISTANCE) 

In this thesis, the technical specifications of the MF Oppedal are adapted to our 
computer simulation, along with the addition of a battery and SC as the HESS. A 
photograph of the MF Oppedal is shown in Figure 8 and its key specifications are listed in 
Table 4. 



Figure 8. MF Oppedal. Source [45]. 


Table 4. Key specifications for MF Oppedal. Source: [46] 


Ferry 

Fength (m) 

L PP (m) 

Breadth (m) 

Draught(m) 

Velocity 

(Kts) 

MF Oppedal 

113.96 

104.9 

16.8 

3.36 
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1. Drag (Resistance) 

Any vessel underway will experience a net force opposing its forward movement 
due to the stress and shear forces acting on the hull surface [23]. This resistance is also 
known as drag. Most of the ship fuel is used to generate the required propulsive power to 
overcome this drag in order to propel forward. As shown in [23], Equation 4 shows the 
general resistance based on vessel design speed (v), wetted surface area (S ), fluid density 
(p ), and the resistance coefficients, (c r ). 


Rt ~ 2 Pseawater V $ C T ( 3 ) 

The total resistance coefficient c T is given by 

c T =(l + k)c F +c R +Ac f +c m (4) 

where k is the form factor, c F is the frictional coefficient of the vessel, c R is the residual 
resistance coefficient (graphical determination according to Harvald and Gulhammer [47]), 
A c F is the roughness allowance, and air resistance coefficient, c AA [23]. 

The form factor is given by 

k = -0.095 + 25.6(c b / ((L w; / B) 2 + (B / T) (1/2) ) (5) 


where L wt is the length of the ship on the waterline, B is the breath of the ship, and T is 
the draught of the ship. As explained in [48], the block coefficient, c B , of the ship is given 
by 


c 


B 


L ,BT 

wl 


( 6 ) 


which depicts the relationship between hull displacement V and the volume in m 3 , 
representing the ship weight. The frictional coefficient according to the ITTC-57 model 
correlation line highlighted in [49] is calculated by 


0.075 

(log 10 Re-2) 2 ’ 


(7) 
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where the Reynolds Number Re is calculated in relation with the viscosity 77 [23] as 

vL , 

Re = —^. ( 8 ) 

n 

In order to determine c R from the Harvald and Gulhammer diagram [47], the length 
displacement ratio (LDR) is calculated in (9) as described in [48]: 

LDR = —(9) 

y/3 

As explained in [50], the Froude Number/^ is a dimensionless parameter measuring the 
ratio of the inertia force on an element of fluid to the weight of the fluid element. It is used 
to find the resistance coefficient c R via the corresponding Harvald chart [47] 


F = 



( 10 ) 


From [22] and [52], the formula for the wetted surface area S is 


5=0.99 


4 V 

- + 0.19 LJ 


( 11 ) 


Equation 12 [23] represents the ratio between volumetric displacement and the 
product of L wl with the cross-sectional area of the submerged hull A m \Al\. 


C " (LA.) 


( 12 ) 


The air resistance coefficient [23] is determined by 


c - c. 


f n A ^ 

Pair A x 


vA 


seawater J 


(13) 


where c x is the air-drag coefficient of the vessel above the waterline, and 0.8 is usually 
used as the default value [49]. A x is the projected area of the ship above the waterline. 

The hydraulic roughness of the hull, also known as the roughness allowance 
coefficient, A c F is determined by [49] 
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+ 0.000125, 


(14) 


A c F 


0.044 



V L wl J 


-lORe 


where k s indicates the roughness of the hull's surface, and the standard value of k s is 
150* 10' 6 [49], 

Equation 15 [23] shows the effective power required to drag the ship through 
the water. 


P e =R t v (15) 

However, we need to consider the efficiencies from the propellers rj 0 , transmission rj s , and 
hull r] H and relative rotative efficiency r] R . Hence, Equation 16 shows the necessary power 
P , based on the efficiencies. 

P P 

P p = — = 7-*-x (16) 

It (VoPsVhPr ) 

For this thesis, the total efficiency r] T for this ship is assumed to be 70%. P is also 

often referred to as mechanical brake power P B [52]. Based on Equations 3 to 16. Figures 

9 and 10 show that as ship speed increases, the drag of the ship increases exponentially as 
along with the mechanical break power required to drive the ship forward. 
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Figure 9. Drag versus speed graph. 
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Figure 10. Brake power versus speed graph. 
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2. Ship Generators 

The MF Oppedal ship is designed to travel up to 13 kts, and has four diesel engines 
aboard [23]: 


• Two Mitsubishi GenSet S12RMPTA with 1100 kW each 

• Two Mitsubishi GenSet S6R2MPTA with 640 kW each 

B. LOAD MODEL 

In our simulations the wave model simulates the load changes experienced by the 
ship sailing under various sea conditions. To simulate the wave action experienced by the 
ship, a sine wave was used to represent ship power demand. The sine wave would change 
in magnitude to simulate the changes in the sea condition from a light varying condition to 
a heavy varying condition as if the ship is sailing through a very severe sea state condition. 

Although in this thesis we used only one stationary frequency component, the 
approach we present is general enough that it can be extended to more complex models 
with a number of frequency components that might even change with time. 

As the sea conditions worsen, more power is required, and thus, a higher magnitude 
of the sine wave. Three wave profiles (15%, 25%, and 45%) are simulated for this paper. 
The first profile is a low-profile load simulating a small sea state. This implies that apart 
from the power P B needed to keep forward motion at a constant velocity, an additional 

15% of power is required to compensate for the effect of wave. The next profile is to 
simulate a moderate sea condition at 25% of the propulsion power and 45% for simulating 
heavy sea conditions. 

The frequency of the profile depends on the wave, usually measured by when the 
two crests of a wave pass a point within a certain time frame. Adopted from a previous 
work [11], in our simulations we assume the wave period is 8 secs, which yields a 
fundamental frequency of 0.125Hz. 
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C. HYBRID ENERGY STORAGE SYSTEM 

Based on the characteristics of the different energy storage devices shown in 
Figure 2, we can determine that no single device type can respond optimally to both high- 
and low-frequency power exchanges. The most common way to address this limitation is 
by combining multiple types of energy storage devices, thus forming the hybrid ESS (or 
HESS). A battery-SC combination is considered in most HESS developments because of 
the availability of these components, their similarity in working principle, their relatively 
low cost, and most importantly, the ability of each component to mitigate the limitations 
of the other component effectively. This combination is also commonly used for 
electrically-driven hybrid vehicles and residential energy storage microgrids and has been 
discussed extensively in much of the scholarly literature [30], [33],, [56]-[60]. 

Before establishing the battery-SC HESS size, we consider the specifications of the 
lithium-ion battery and supercapacitor, as shown in Tables 5 and 6. 


Table 5. Specifications for lithium-ion battery. Source: [59]. 


Lithium-ion Battery 


Nominal Voltage 

173V 

Nominal Capacity 

260AH 

Energy 

45kWh 


Table 6. Specifications for SC. Source: [60]. 


125V Supercapacitor 


Rated Voltage 

125V 

Rated Capacitance 

63F 

Maximum Continuous Current 

140 A 

Usable Specific Power 

1,700 W [60] 


The battery-SC HESS is sized according to the load generated by the highest wave 
profile of 45%. There can be many permutations for the HESS. The arrangement of the 
battery bank is designed to meet the power system requirements by adding batteries and 
SCs in a series to reach the desired direct current (DC) voltage. 
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D. HYBRID POWERTRAIN ARCHITECTURE 


The powertrain of the hybrid ship considered in this thesis incorporates generators, 
batteries, and SCs. The total power required by the ship is split between the generators and 
the HESS. The generator will provide the steady current based on the desired speed, and 
the HESS will take up the varying current required due to the wave load profiles. 

The power control presented in this thesis is based on a simplified dynamic model 
of the ship moving in the forward direction, propelled by its own power and subject to 
external forces such as drag and wave action. Its dynamics then can be written as 

Mv(t) = F pwp (t)-F drag (t)-F wave (t), (17) 

where F is the propulsion force, F, is the drag force, and F wave is the force due to the 
wave. In steady state, the forward velocity is expressed as 

v(0 = v+v(0 , (18) 


where v is the steady state forward velocity and v(t) is the velocity perturbation due to the 
wave. Then, expanding the drag resistance in a Taylor series around the steady state 
velocity v , we obtain 

mt) = F prop (0 + F prop (t) - F drag (v ) - VF dmg (v )v(t) - F wave (t) (19) 


The following equation 


prop 


= F V 

prop 


( 20 ) 


is the power required to move the ship forward at the steady speed v . 

For the wave component, the propulsion power required to compensate for the wave 
action and keep the velocity perturbation v(r) to a minimum is 


P =F (t)v 

prop prop v / 

x F prop {t)v 

Combining Equations 19-21, we obtain 


(21) 
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(22) 


V(t) _ .c v(t) , 1 

+ ^-(P pr o P (t ) ~ P wave ( f )) ’ 

v v 2 (a) 

where S - VE ira? (v) / M and /f = — MV 2 , where M is the mass of the ship. 

Since the propulsion power comes from a number of energy sources, we can 
express that 

P prop (t) = rj eff (P gen (t) + P batt (0 + 4 ( 0 ) (23) 

with ri eff as the efficiency of energy conversion from electric to mechanical energy. 
Assuming a constant bus voltage, we yield 

v(t) v(t) V ~ 

= -S^ + ^% ff (I gm {t) + I batt {t) + l sc {t)-I wav M ( 24 ) 

v v 2(Z) 

with I wave (t) denoting the current required to fully compensate for the wave action on the 
forward velocity. Define 


x v (0 = 


2 (K) v(t) 

^bus Veff V 


and Equation 24 becomes 

x v (t) = -8x v (t) + I gen (t) + I batr (t) + I sc (t) ~ T wave (t). (26) 

Expressing this equation in discrete time by approximating the derivative, we 


obtain 


(t + 1 ) = a ■ - ST )X V ( t) + l gen (!) + I battery (0 + I sc o t ) - I wave (f ), 


where t is now the discrete time index and 


x v (t) = 


x v (t) _ 2 (K) m 
Ts V bus T s v 
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In these equations the bus voltage of the ship is 800V and assumed constant. If we 
define I B as the steady state current to move the ship forward at a constant velocity v the 

total current required is the sum of the wave current I wave and the brake current I B . 

1. Diesel Generator Model 

The diesel generator is also modeled in discrete time. The generator has two current 
components, I B and I . I B provides a steady-state current based on the required power 

P B necessary to keep the ship moving forward and maintaining its course. For this thesis, 
the goal of the MPC controller is to provide current, I , for efficient use of the generator 

and energy storage sources in the presence of adverse sea conditions. Since it is desirable 
to keep the generator operating at a steady pace, avoiding deep swings in power, we define 
the state space model of the generator in terms of the current delivered I , where A/ 
is the control variable, constrained to small rates of change: 

I g en(t + V = I g en(t) + U gen (.t). (28) 

2. HESS Model 

The following section will elaborate on the state space equations for the battery and 
SC. 


a. Battery 

The main parameters of the battery are its capacity in Amps per hour (Ah) and State 
of Charge (SOC). The rated capacity is simply the energy capacity of the battery under 
normal conditions. For example, 200 Ah means if the battery is fully charged, in principle 
it could provide a sustained current of 200 Amps for one hour before it is completely 
discharged. The SOC reflects the level of charge of the battery relative to its capacity. The 
unit of SOC(t) is a percentage, and it indicates the amount of energy stored at time “t” 
versus the energy stored at rated capacity. There are also additional parameters defining 
battery performance such as the State of Health, State of Function, and shelf life. Those 
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factors will not be considered for this thesis. Mathematically, the battery state-space 
equations are defined in Equations 29 and 30. 


I battery it ^ 1 I battery (o+awo ( 29 ) 

soc Batt (t+ 1 ) = soc Ball (0 - rhatt (0 (30) 


Adapted from [31], y represents the charging and discharging factor of the battery. 
A positive current represents that the battery is providing current. The Coulomb counting 
method is used to determine y, which can be defined as 


r = 


3600 (Q ) 


(31) 


b. Supercapacitor 

The fundamental properties for the SC are its voltage and capacitance, as the energy 
of a capacitor is defined as [61]: 

E = -CV 2 (32) 

2 

Similar to the battery, the SOC(t) of the SC has the same meaning as energy stored 
at time t versus energy stored at rated conditions. Mathematically, the state-space equation 
for the SC current and SOC are shown in Equations 33 and 34. 

I sc (t + \) = I sc (t) + AI SC ( t ) (33) 

SOC sc (t +1) = SOC sc it) ~ sc (0 > (34) 

where the [3 represents the charging and discharging factor for the SC. Adapted from [31], 
the value /3 can be defined as 


P = 


C V 


(35) 


c. State Space Model and Constraints 

The state variables and inputs of the optimization problem are defined by 
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x(t) = 

X v 9 1 gen 9 1battery 9 SOCba tt ery > I SC ’ ^OC sc 


~lT 


(36) 


and 


n(0 = [ A/ ?eil , A I batt , A/J. 


(37) 


In addition, we define the output of the system as y(t) = I residua , (0 , which represents 

the difference between the total current generated and the current required to fully 
compensate for the wave actions. This can be written as 

I residual i f ) = ^gen i f ) + halt i f ) + he (0 _ have (0 ■ (38) 

Equations 27-30 and Equations 33-34 can be consolidated into a state-space matrix 
in the following form 


x(t + 1) = F.x(t) + Gu(t ) + Bw(t ) 


1 residual^ = #*(0 + MO + Dw i f ) . 


(39) 

(40) 


where 


F = 


"1 1 

1 

0 

1 

0 " 


“0 

0 

0 " 


'-1 

0 “ 

0 1 

0 

0 

0 

0 


1 

0 

0 


0 

0 

0 0 

1 

0 

0 

0 

G = 

0 

1 

0 

B = 

0 

0 

0 0 

-y 

1 

0 

0 


0 

0 

0 


0 

0 

0 0 

0 

0 

1 

0 


0 

0 

1 


0 

0 

0 0 

0 

0 

-fi 

1 


0 

0 

0 


0 

0 

= [o 

1 1 

0 

1 

0 ] 

J = 

[0 

0 

0 ] 

F> = 

[-1 

0 ] 


At any time t, assuming we have full state observation x(t), the goal of the 
controller is to find a sequence of predicted future inputs u lU ,u lMk ,...,u l+NU , which 
minimizes the cost function 


, 1 N-\ 

V (M),w +tt ) = +u t + ku T Ru t + k\,) + { x m T Eu t\k) 

7 ^ k=Q 


’ ( 41 ) 
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where x t+a indicates the predicted state at time t + k > t and x At = x(t) the current state. In 

this expression the rightmost term expresses the total cost from time t + N to infinity, and 
the corresponding matrix P is computed from the Riccati equation associated with the 
Discrete Linear Quadratic Regulator. 


In the application of this thesis, the parameters associated with the cost matrices Q 
and R are expressed in terms of the respective state variables as 

v ( x «’ H „„) = yZ(<7 x„ I V lf+M Q Igen |-^ ^ |-^ Q SOCbatt \sOC bam \ M + qJlscL + qsoJSOCscl 


1 

+ — 
2 


X(r. V A/ ) (.v(/) /:.///(/) J-v/;,,/ J v ( . Ai( . (42) 


The ’r’ and 'q 1 subscripts denote the weighting factors on the states and input. Table 
7 highlights each subscript and its function. 
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Table 7. Weightage factors on the states and input. 


Subscript 

Descriptions 

q 

-L XV 

Weightage on the x v 

The weightage determines 
how much emphasis is 
placed on the energy 
devices states; i.e., a higher 
factor weightage equates to 
a high penalty on the state. 

E g-, if <7 is » than n , 

1 Ibatt 1 Isc 

the SC will be utilized more 
than the battery. 

^ Igen 

Weightage on the j 

4 Ibcitt 

Weightage on the j hatt 

G SOCbalt 

Weightage on the Battery SOC 

q ISC 

Weightage on the J s( . 

G SOCsc 

Weightage on the SC SOC 

V gen 

Weightage on the rate of change 
for A T 

J- gen 

This controls the rate of 
change for inputs. 

V batt 

Weightage on the rate of change 

fOT A Ibar t 

Vsc 

Weightage on the rate of change 

for AT 

■L SC 


In order to better illustrate the specific implementation of the MPC, we follow an 
example where the time horizon has a length N=3 samples. Again, defining t+klt =t+k, with 
k>=0 integer, the predicted time instant given the present time t, Equation 43 yields a 
simple recursion 

A+u, = Fx(t) + Gii(t) (43) 

X t + 2\, = Fx , +a + Gu r +U = F2x (0 + FGu t+ l\, + Gu (0 ( 44 ) 

X t+3\t = Fx t+2\t + Gu , + 2 = E 3 *(0 + F 2<Gu , + 21, + FGu t+ w + Gu (0 • (45) 
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In this recursion, x(t) indicates current state, which is observed, and x t+klt ,k > 0 
indicates predicted states, while u t+k]t for k>=0 denotes values of the inputs to be 

determined by optimization. In the MPC framework only u(t)t = u m will be the actual 
input to the system. 

The three preceding equations can be rewritten in block matrix form as 


x(t) 


7 


0 

0 

0 " 

u(t ) 

x t+\\t 


F 


G 

0 

0 


= 


x(t) + 

FG 


0 

M r+ii/ 

X t+2\t 


F 2 

G 

X t+3\t 


F 2 


F 2 G 

FG 

B 

U t+2\t _ 


= S x x{t) + S u u t 


The terms in the cost function introduced in Equation 43 can then be written as a 
standard quadratic form as follows: 


Qx t+k „) + (x t+3U T Px , +3{ ,) = [ x(t) T 

k=0 



Q 

0 

0 

0' 

~m T 

r 1 

0 

Q 

0 

0 

x t+l\t(f) 

X t+3\t J 

0 

0 

Q 

0 

X ,+2\t (0 


0 

0 

0 

p 

X l+3\t (0 


=x T (t)S T x QS x x(t ) + u T t S T U QS u u t + 2x t ( t)S T x QS u u, 


while the other term of the cost function can be written as 


JV-l 

'^j{ U t+k\k ^ U t+k\k 
k =0 


)-[ 


' U t\t U t + \\t 


*f+2lf 


R 0 0 
0 R 0 
0 0 


t+i\t 


t+2\t 


= u , Ru, 


Putting everything together yields 

V ( x(t),u t ) = ujY,u t + 2 x(t) T Om, + x(t) T 'Vxit ), (46) 
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where r = S' u QS H + R + S 7 E + E T S u , ® = S T X QS U + S u E and T = S T OS x .This function is 
optimized with respect to the input vector u t . The control input at time t then is the first 
element of the vector u t , as u(t) = u rit . 

d. Setup Constraints 

The goal is to set up the minimization of the cost function V(x(t), u(t)) subject to 
the constraints that are based on the design consideration of the energy storage devices. 
Hence 


-20A < X v < 20A 

(47) 

-200A < I batt < 200A 

(48) 

20% < SoC bnn < 90% 

(49) 

-240A < I sc < 240A 

(50) 

10% < SoC sc < 90% . 

(51) 

Compiling Equations 47 to 51 into the matrix form, we get 

A x x(k)<b x k=l,...N 

(52) 

A Ak)<b u k = 1,...N. 

(53) 


For example, the constraints just defined for the battery current and SOC could be 
written in matrix form as 

"1 0 ] [ 200 " 

-1 0 200 

x(t) < 

0 1 0.9 

0 -lj |_-0.2 

The aim of the optimization problem is to minimize the excursion of the generator 
current I gen while compensating for wave action, subject to the constraints of the energy 

storage elements in terms of their currents I ball , I sc , and their states of charge. To this end, 
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we make use of the MATLAB function mpcqpsolvers, which solves a quadratic 
programming problem in the form of Equation (42) with constraints using the KWIK 
algorithm [62]. 


e. Francis Equation 

Since the model includes an external signal w(t ), which models the effects of wave 
on ship power, a particularly suitable approach is given by the use of the Francis Equation 
introduced in [63]. This approach is based on the fact that wave can be assumed as periodic, 
at least within a finite time interval, represented by sums of sinusoids. Consequently, they 
can be expressed in a state-space model which, in the case of one sinusoidal component, is 
of the form 

w(t + \) = A w w{t), (54) 

where 


cos 0 -sin 0 
sin 0 cos 0 


(55) 


and 9 being the digital frequency of the wave and the state w(t) = [vq, vv 2 ] 7 . The wave 

act as an additional force that the ship has to overcome in addition to the drag, which was 
highlighted in Equation 39. It can be shown that steady state inputs and states solutions of 
the form 

u(t) = Lw(t) (56) 


and 


x(t) = Twit) 


(57) 


can be computed with L and T matrices of suitable dimensions (Figure 11). 
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Model 


x(t) = Tw(t) 

> 

► 


u(t) = Lw(t) 


Figure 11. Illustration of the Francis equation 


The matrices T,L are solutions of the Francis equation defined as follows: 


FT + GL + B = TA, 

(58) 

HT + JL + D - 0 

(59) 

It is shown in [63] that a solution exists provided the matrices A and F have no 

common eigenvalues. Furthermore, based on the preceding definitions, we define new 

states z(t ) and input v(f), respectively, as 

z(t ) = x(t)-Tw(t ), 

(60) 

II 

1 

Fh 

(61) 

Substituting Equations 61 and 62 into Equations 1 and 2, we obtain 


z(t + 1) = Fz(t) + v(t) 

(62) 

and 

I r Jt) = H Z (t) + Jv(t). 

(63) 


The goal is to find a control input v(t) for regulation of the system, so that the 
corresponding input u (/) = v(f) + Lw(t) regulates the system. 


y(t) 
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Since the constraints stated in Equations 47 to 51 refer to the state and input vectors 
x(7)and //(/), we need to determine constraints in terms of z(t) and v(t). Recalling the 

constraints in x(t )and u(t) domain 

A x x(t)<b x (64) 

and 

An(t) < b u , (65) 

substituting Equations 64 and 65 into Equations 52 and 53, which yields 

A x z(t) < b x - A x Tw(t) (66) 

and 

A u v(t) < bu - A u Lw(t ). (67) 
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IV. RESULTS AND ANALYSIS 


This chapter presents the results obtained after running several time-domain 
simulations of the model. No internal losses are assumed for the dynamics and they are 
assumed to be operating under ideal conditions. The aim of the model is to present the 
MPC algorithm and how it responds using different values of the weights in the cost 
function defined in Table 7. The model is implemented in MATLAB for three wave 
profiles. In particular, each scenario also compares the battery and SC current outputs based 
on different values of the weight q Ibatl in Table 7. 

The MPC model illustrated in Chapter III is evaluated against six values of the 
weight q Watt (1, 5, 10, 20, 30, and 50). The objective is to keep the generator running at a 

steady state while the HESS is compensating for the wave when the ship is travelling at 
llkts. The constraints stated in Equations 47-51 are also applied in the model. The 
simulation time for the model is 28,800 seconds, or 8 hours, so as to ensure the model is 
able to reach a steady-state output. 

A. LIGHT SEA CONDITIONS (15% WAVE PROFILE) 

In light sea conditions, the additional power required to navigate through the wave 
is approximately 390 kW, which is 15% of the power required to cruise in calm seas. The 
MPC model is able to optimize and minimize the generator current, / , to a steady state 

of 108 A for all q Ibalt , while the battery and the SC are providing the rest of the required 
current to compensate for the wave action and keep the forward velocity of the ship as 
smooth as possible. Figure 12 shows the residual current, I residuai , at steady state and the 
expanded view of the initial transient response. 


35 





0 20 40 60 80 100 

Time(sec) 


(a) 


(b) 


Figure 12. (a) I residual with q Ibatt and q Isc = 1. (b) Enhanced view of the I residual 

initial transient response. 


Figure 13 shows the velocity perturbation v(t ). The model can reach steady state 
quickly and the velocity fluctuation is kept small. The goal is to keep the residual current 
I residual minimal. It is not possible, however, to minimize v(t) and I residual to zero due to the 
constraints. 
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Figure 13. (a) v(t) with q Ibatt and q Isc = 1. (b) Enhanced view of the initial 

transient response. 


As already observed, Figure 12 shows the plot for I residual , where it can be noted that 
the transient response is small compared to the two other bigger wave profiles. The current 
I residual oscillates between +13A at steady state for the rest of the simulation time. 

Figure 14 shows the currents plots for I wave , I , I hatt , and I sc . Figure 14(a) shows 

that the total current generated by I , I hatt , and I 5C is able to compensate for I wave and 
with the slight difference is as shown in Figure 12, while Figurel4(b) shows the individual 
current of I gen , I batt , and I sc 
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Figure 14. (a) Enhanced view of the I gen , I hatt , and I sc output currents with 

q Ibatt and q lsc =1. (b) A dissected view of the individual current. 


As the value of weights q Ibatt and q kc are the same, the output currents from both 
devices are approximately the same. The slight difference is due to the different y and P 
factors related to their SOCs. The rate of change for the generator A I was set to 50 A 
per sample, corresponding to 62.5 Amps/sec to illustrate the slower rate of change as 
compared to the higher A I batt , and A I sc , and the difference is shown in Figures 14 and 16. 

As the I wave is small, the rate of change for I is not that significant and this observation 
can be seen in the higher wave profiles where the I wave is larger. It is also observed that 
/ gen is slightly slower than the other two currents due to A/ . 
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Figure 15. 


(a) 
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residual with 


y ,hatt ^ and ^ Isc '. (b) Enhanced view 
initial transient response 


I 


residual 


Figure 15 shows the I residua i when q Ibalt is increased to 30. It is observed that the 
transient response has a higher peak value although it has approximately the same settling 
time as Figure 12 The I residua i can be further reduced by limiting the range of x v . For this 

thesis, x v is unchanged in order to observe the impact of q Ibatt . I residual is able to reach 
steady-state and maintain at +16A throughout the simulation time. 
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Figure 16. (a) Enhanced view of the I gen , I batt and I sc output currents with 
Qibatt = 30 and q hc = 1. (b) A dissected view of the individual current 


Figure 16 shows the current plots for I wave ,/ , I batt and I sc . Figure 16(a) shows 

that the total current generated by I , I batt , and I sc is able to compensate I wave with the 
slight difference as shown in Figure 15, while Figurel6(b) shows the individual current of 
I S m > han > an( 3 I sc . As the values of the weights increase, the difference in current between 
the battery and the SC increases, and as compiled in Table 8, the current differences start 
to taper off when q Ibatt is greater than 20. 

Table 8 compares the values of the I gen , I ban , and I sc peak output currents for the 
six different values of q Ibatt . As the I wave is considered small, the difference between each 
q Ibalt is small. The difference would be more apparent as the wave profile increases. 
Nevertheless, it is observed that the I batt decreases as q Ibatt increases. It is also observed 
that when q Ibatl is greater than 20, the difference between the battery and the SC starts to 
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taper off. Similar observations are seen for the larger wave profile. It is also noted that I 
remained the same for all values for q Ibatt . 


Table 8. Compilation of the peak output currents for 1 , I hau , and I sc . 


Qlbatt 

I gen 

batt 

L 

1 

107.9A 

188.63A 

188.83A 

5 

107.9A 

175.17A 

200.8A 

10 

107.9A 

171.23A 

204.47A 

20 

107.9A 

168.76A 

206.79A 

30 

107.9A 

167.84A 

207.66A 

50 

107.9A 

167.06A 

208.39 


B. MODERATE SEA CONDITIONS (25 % WAVE PROFILE) 

The moderate sea condition is commonly encountered by a ship during deployment, 
and wave amplitudes for this condition are slightly larger as compared to the light sea 
condition. Therefore, the ship is required to generate more power to compensate for the 
wave in order to maintain its velocity. The ship requires an additional 650 kW to navigate 
through the wave. For all values of q Ibatt , the MPC controller is able to optimize the states 

and minimize the generator peak output current I to 163 A, while the battery and SC 
are providing the required currents to compensate for the wave. Similar to the light sea 
condition, v(r) is observed to be small and this is consistent for all values of q Ibatt . 

Figure 17 shows the velocity perturbation v(t) during the moderate sea condition. 
As compared to the light sea condition, the model took a longer time to reach steady state, 
while the fluctuation is small. A similar observation is seen for all values of qbatt. 
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Figure 17. (a) v(t) with q Ibatl and q hc = 1. (b) Enhanced view of the initial 

transient response. 


Figure 18 shows the current difference I residual between the I wave , I gen , I batt , and 
I sc . Similar to Figure 17, the model took slightly longer to reach steady state as compared 
to light sea conditions. The steady state I residua , is also slightly larger as compared to light 
sea conditions. This is expected as the required I wave is larger. The I residual stabilized 
between +18A for the rest of the simulation time. 

Figure 19 shows the currents plots for I wave , I , I batt , and I sc . Figure 19(a) shows 
that the total current generated by I , I batt , and I sc is able to compensate for I wave with 
the slight current difference is as shown in Figure 18. Figure 19(b) shows the individual 
current of I gen , l hau , and I sc . 
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Figure 18. 


(a) I residual with Qoatt and c Usc = 1 • ( b ) Enhanced view of I residual 


initial transient response. 
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Figure 19. (a) Enhanced view of the / , I batt , and I sc output currents with 

q Ibatt and q hc = 1. (b) A dissected view of the individual current. 
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As the value of the weights q Ihatt and q hc are the same, the output currents from 
both devices are approximately the same. The slight difference is due to the different y and 
(3 factors. The rate of change for the generator A I remains unchanged from the previous 

sea profile so as to illustrate the difference in rate of change as compared to the battery and 
the SC. As compared to light sea conditions, this is more evident in Figures 19 and 21 as 
/ waveform is observed to be slower than the rest of the current waveforms. 



Figure 20. (a) I residual with q Ibatt = 30 and q hc = 1. (b) Enhanced view of the 

/ ra7dHfl , initial transient response. 


Figure 20 shows the plot of the l residuai when q Ihatt is increased to 30. It is observed 
that the transient response is slightly longer than that shown in Figure 18 when q Ibalt is 1. 
This is because a higher penalty is put on the battery, thus reducing the battery output 
current at each sampling time. The I residual is able to maintain the steady-state oscillations 
between +22A throughout the simulation time. 
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Figure 21 shows the currents plots for I wave , I gen , I hatt , and I sc when q Mlt is 
increased to 30. Figure 21(a) shows that the total current generated by I gen , I hatt , and I sc 
is able to compensate I wave with the slight current difference shown in Figure 20. Figure 
21(b) shows the individual current of I gen , I batt , and I sc . 



Figure 21. (a) Enhanced view of the/ , I hatt , and I sc output currents with 

c h>, u n = 30 and q Isc = 1. 


As the value of the weight increases, the difference between the currents delivered 
by the battery and the SC, respectively, increases. Similar to the previous wave profile of 
15%, the current differences start to taper off when q Ibatt is greater than 20. In comparison 

to the light wave conditions, the current differences for the q Ibatt in moderate sea conditions 
are larger. Table 9 compiles the steady-state peak output current for I , I batt , and I sc 
based on different values of q Ibatt . 
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Table 9. Compilation of the peak output current for 1 , I hatt , and I sc . 


Q-Ibatt 

I gen 

1batt 

Isc 

1 

163.17A 

344.49A 

345.01A 

5 

163.16A 

305.60A 

385.77A 

10 

163.15A 

293.53A 

397.98A 

20 

163.14A 

285.78A 

405.67A 

30 

163.14A 

282.85A 

408.54A 

50 

163.13A 

280.36A 

410.96A 


In this simulation, I remains at 163 A while the battery and SC are providing 
the required current to compensate for the effects of the wave expressed by I wave . The 
difference in current between the battery and the SC is more substantial as q Ibatt increases. 
With the slower rate of change of the generator, the battery and SC is able to compensate 
the I wave response effectively. This reduces the loading stress on the generator and reduces 
defective occurrences due to loading stress. 

C. HEAVY SEA CONDITIONS (45 % WAVE PROFILE) 

The heavy sea condition is rarely encountered, and it represents a ship navigating 
through severe weather at sea. Under such a situation, without the ESDs, the generator will 
be placed under immense stress from generating the required fluctuating power when 
navigating up and down the wave. For this condition, the ship requires an additional 1.17 
MW of power in order to compensate for the wave action and keep the forward velocity as 
steady as possible. As the demand is much higher, the rate of change A I for this 

simulation also increases to 250 Amps/sec, indicating that more generators are operating 
to provide more power as well as to spread out the loading stress. After optimizing the 
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states and inputs by the MPC controller, the generator output current, I gen , is maintained 
at 398 A, while the battery and the SC are generating the required currents to compensate 
for the wave. Similar to the previous two wave profiles, different q lban weightages are 
analyzed. 

Figure 22 shows the velocity perturbation v(r) during the heavy sea condition. As 
compared to the previous two sea conditions, the model has a bigger initial transient 
response and v(r) is larger. With the increased value of A/ the model is able to reach 

steady state quickly. A similar observation can be made for all values of q Ibatt . 
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Figure 22. (a) v(t) with q lbatt and q hc = 1. (b) Enhanced view of the initial 

transient response. 


Figure 23 shows the current I resjdual , which is again the difference between the I wave 

and the total generated current, which is the sum of I , I hatt , and I sc . The transient 

response has a higher peak value than the previous cases. The model is able to reach steady 
state quickly and is maintained at ± 30A for the rest of the simulation time. Increasing the 
A/ will reduce the transient response and the steady state difference, but that would also 

indicate either more generators are connected to the ship microgrid or the generators are 
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taking on additional loading stress. An alternative is to select a battery and SC that have 
much larger capacity than the benchmark. 
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Figure 23. (a) 


residual with ^ ffia "and 

transient response 


^ Isc 1. (b) Enhanced view of ^ residual initial 


Figure 24 shows the currents plots for I wave , I , I batt , and I sc . Figure 24(a) shows 
that the total current generated by I , I batt , and I sc is able to compensate for I wave with 
the current differences as shown in Figure 23. Figure 24(b) shows the individual current of 
/ , I han , and I sc .As the values of the weights increase, the difference between the battery 

and the SC increases as well. Similar to the previous two wave profiles, the current 
differences start to taper off when q Ibatt is greater than 20. 
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Figure 24. (a) Enhanced view of the I gen , I batt , and I sc output currents with 

q Ibalt and q Isc = 1. (b) A dissected view of the individual current. 


Figure 25 shows the I residual when q Ibalt is increased to 30. It is observed that the 
transient response is much larger when q Ibatt is 1, as shown in Figure 23. The I residuaI is able 
to reach steady state and maintain oscillations of +42A throughout the simulation time. 
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Figure 25. (a) I midml with q Ihatt = 30 and q Isc = 1. (b) Enhanced view of I residual 

initial transient response 


Figure 26 shows the currents plots for I wave , I gen , l batt , and I sc when q Iball is 

increased to 30. Similar to the previous two wave profiles, as the weightages increase the 
difference between the battery and the SC increases. The current difference between the 
battery and the SC, respectively, starts to taper when q Ibatt is greater than 20, while the / 

is maintained at 398A. With a higher A I , the ESDs are still able to generate most of the 
required current in order to compensate for the wave current. 
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Figure 26. (a) I gen , l baU and l sc output currents with q Ibatt = 30 and q lsc = 1. 

(b) A dissected view of the individual current 


Table 10 shows the steady-state peak output current for the I gen , I han , and I sc 
based on different q Ihatt weightages. The difference between the currents is not as 
substantial it was for the moderate sea condition due to the increased A/ . 


Table 10. Compilation of the peak output current for / , I batt , and I sc . 


Q-lbatt 

I gen 

^bait 

he 

1 

398.20A 

532.48A 

532.79A 

5 

398.46A 

511.89A 

546.47A 

10 

398.50A 

505.91A 

550.18A 

20 

398.52A 

502.13A 

553.13A 

30 

398.53A 

500.70A 

554.31A 

50 

398.52A 

499.45A 

555.32A 
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One common observation for all three wave profiles is that the battery and the SC 
output currents begin to taper off when q Ibatt is greater than 20. This highlights that an 
increase in the difference between the states greater than 20 does not have a significant 
impact for this model. The ESDs are able to compensate for most of the I wave , thus reducing 
the stress on the generator significantly. 

D. ESD STATE OF CHARGE 

Figures 27 and Figures 29 compare the battery and the SC SoC at two different 
values of the weight q Ibatt ( q lbalt = 1 and q lbatt = 30 ), respectively, using the moderate wave 
profile as this is the scenario most commonly encountered during sailing. It is observed 
that the SC SOC reaches its maximum rated value of 90% more quickly at a lower q lbatt as 

compared to a higher q Ibatt weightage. This is expected as both ESDs are providing 

approximately the same amount of current, and when the SC has reached its rated capacity, 
the additional current during the charging phase is diverted to charge the battery. This is 
evident as the battery SOC starts to increase when the SC SOC has reached its upper limit 
set in Equation (51). As q Ibatl increases, the SC is providing more current as compared to 

the battery. The SC SOC requires more time before reaching its maximum rated SOC and 
before the additional current is diverted to charge the battery. 
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Figure 27. (a) Battery SoC. (b) SC’s SoC when q, hatt = 1. 
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Figure 28. Enhanced views of (a) Battery SoC (b) SC SoC when q Ihatl = 1. 
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Figure 29. (a) Battery SoC (b) SC SoC when q Ihatt = 30. 
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Figure 30. Enhanced views of (a) Battery SOC (b) SC SOC when q Ibatt = 30. 
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V. CONCLUSION AND RECOMMENDATIONS FOR FUTURE 

WORK 

This chapter presents the conclusion based on the results of the model developed 
and tested in this thesis. Areas for future work are also recommended, including 
implementing the model for different kinds of ships, designing a Kalman filter to estimate 
the state of the dynamic model for the wave action, having different charging and 
discharging factors for the battery and SC, and simulating a more realistic wave model. 

A. CONCLUSION 

In this thesis, we presented a control system that uses a combination of power 
generation and energy storage to keep a ship on a steady pace in different sea conditions. 
Although the work is preliminary and based on ideal, simplified conditions, it shows that 
a control system that uses the Francis Equation [63] to model periodic perturbations, and 
Model Predictive Control (MPC) to include physical limitations, can be a viable solution 
to this problem. 

In particular, it shows that the integration of ESDs in the ship microgrid provides 
numerous advantages, such as increasing fuel efficiency, reducing loading stress on the 
generator, ensuring stability in the overall powertrain, and providing cost benefits due to 
improved generator maintenance. Moreover, it also helps provide environmental benefits 
as the generators are operating at steady-state speed with the aid of the ESDs. 

A model of this sort could also be developed as a useful tool that enables the 
designer to size ESDs, such as batteries and supercapacitors, for practical applications. This 
conclusion is based on the fact that, at each iteration, the MPC controller searches for an 
optimal solution based on physical constraints of the various devices. If the components 
are not properly sized (for instance, inadequate storage devices) the controller fails to find 
a solution, which is remedied by properly increasing storage capacities. 
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B. FUTURE WORK 


This research serves as a preliminary study into how MPC is an effective energy 
control management tool for a ship microgrid and further expansions can be implemented 
to improve its robustness. 

There are several areas of follow-on research that can be pursued from this thesis. 
One area to explore is to design a Kalman filter to track a varying load profile. For example, 
the ship is navigating at low speed when exiting and entering ports and increasing to 
cruising speed through the open sea. Also, other electrical loads could be considered in the 
simulation. Such simulations might include higher power machineries such as bow 
thrusters operating during docking and weapons systems such as the U.S. Navy 
electromagnetic rail gun. Further, the model could be expanded to other sizes of ships. 

Additional studies could explore different charging and discharge gammas and 
betas for the battery and supercapacitor. Such investigations would improve the model in 
modeling the characteristics of the devices. 

Similarly, the wave model could be expanded to be more realistic. The current wave 
model only represents two-dimensional wave, and that could be extended to more complex 
models with a number of frequency components with might even change with time. 

Lastly, this preliminary study of the MPC algorithm on a HESS could be expanded 
to model a land-based microgrid or forward-operating grid with more than one energy 
source. 
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APPENDIX. MATLAB SCRIPT 


A. MAIN SCRIPT 


% MF Oppedal 

o 

o 


v_kn = 11 % Enter desired speed 

[Pe,v,DSP] = shippower(v_kn); 

v % speed in mete%rs per seconds 
DSP % mass of ship, kg 

Pp = Pe *4 % 4 Engines total requirement 
Vbus = 800; 

I_bus = Pp/Vbus % total current required move 
forward without the current effect, FB-Fdrag 
%% MPC Control 

% Ship Generator Parameters 
Pgwer = 2*(l.le6) ; % MW 

IGeny = Pgwer / Vbus 

% Battery 
Voc =12' 

0 batt = 200 ; % AH, capacity of battery; 
%SuperCapacitor 

% 125V moudle_datasheet_maxwell 125V heavy 

transportation module SC 

Cap_number =7; % 800V/136 = 6 cap 

Vmax = 135; % SuperCapacitor Vmax , Vrated = 

125V 

SC_Cap = 63; % SC capacitance 
%wave parameters 


57 


F_wave = 0.125 ; %Hz 

T_wave = 1 / F_wave ; % 8 Secs 

FS = F_wave*10 % sampling frequency % 

Ts = 1/FS 

random_number = rand*randi([0 2]) * pi; % 

alpha =0%random_number %0.5204 
thetaO = 2*pi*(F_wave/FS) 

%%Setting simulation time 
T_final = 1* 3600 % 1 hour 

N_final= T_final*FS ; 

t = 0:N_final-l; % setting the intergers based 
on simulation time 


% %Variation in load from load profile % 

SS_factor = 0.25 ;% 15% 25% and 45% 

I_wave = (SS_factor*I_bus) * (cos(thetaO*t + 

alpha)); %Iwave based on load profile %Change 
term 

I_Overall = I_bus+I_wave; 


max(I_wave); 
min(I_wave); 

Amplitude = 0.5* (max(I_wave) - 

min(I_wave)) %load_var; % amplitude of the wave 

curve 

%Getting amplitude of the waves , wl and w2 

wl = Amplitude * cos(alpha) % Iwave based on wl 
w2 = Amplitude * sin (alpha) % Iwave based on w2 
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[e, x, u, w]=Powerscript(I_wave,wl,w2,thetaO,FS,IG 
eny,t,Voc, 0 batt,Vmax,SC_Cap,T_wave); 

%% 

% V_tilde_dot / V_bar 

e f f = 0.9; 

K= 1/2 *DSP*(v A 2) 

V_tilde_dot_V_bar = (Vbus / (2*K) )* Ts *eff 

*x (1, : ) ; 

V_tilde = ((Vbus / (2*K) )* Ts *eff *x(l,:) * v 

) ; 


figure (4 9) 
plot(V_tilde) 
title ("V tilde") 
xlabel( "Time (sec)") 
ylabel (" Speed (m/s)") 
xlim([0 length(t)]) 

figure(40) 
plot(V_tilde) 
title( "V tilde") 
xlabel ("Time (sec) ") 
ylabel (" Speed (m/s)") 
xlim([0 100]) 

%% 

set(groot, 'defaultAxesXMinorGrid' , 'on', 'default 
AxesXMinorGridMode' , 'manual' ); 

set(groot, 'defaultAxesYMinorGrid' , 'on', 'default 
AxesYMinorGridMode' , 'manual' ); 

XV = x (1, : ) ; 

Igen = x (2, : ) ; 

Ibatt= x (3, : ); 
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% 


use for compiling 


SOC_Batt = x(4,:); 

ISC= x(5, : ) ; 

SOC_UC = x(6, : ) ; 

x_T = transpose(x); 

%%U(t) 

Delta_Igen = u(l,:); 

Delta_Ibatt = u(2,:); 

Delta_ISC = u(3,:); 

%outputvy(t) 

figure(1) 
plot (e); % y (t) 

title (' Residual Current, I_{residual}' ) 
xlabel ("Time (sec) ") 
ylabel ("Current (A) ") 
xlim([0 length(t)]) 

figure(41) 
plot (e); % y (t) 

title([ 'Residual Current,I_{residual}' ]) 
xlabel ("Time (sec) ") 
ylabel ("Current (A) ") 
xlim([0 100]) 

% Battery Output graphs 

figure (2) 
subplot (2,1,1) 
plot(Ibatt) 

title ("Battery Current (A)") 
xlabel ("Time (sec)") 
ylabel ("Current (A)") 

subplot (2,1,2) 


plot(SOC_Batt) 
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title ("Battery State of Charge") 
xlabel("Time (sec)") 
ylabel ("SOC(%)") 

figure (50) 

% % invert the graph 

SOC_Batt_Inv = SOC_Batt*-l; 
plot(SOC_Batt_Inv) 

title ("Battery State of Charge") 
xlabel("Time (sec)") 
ylabel ("SOC(%)") 


yticks([ -0.08 -0.06 -0.04 -0.02 0]) 
yticklabels ({'72', '74','76','78', '80'}) % set 
the label to percentage 

% % invert the graph , enahnce view 

SOC_Batt_Inv = SOC_Batt*-l; 
plot(SOC_Batt_Inv) 

title ("Battery State of Charge") 
xlabel("Time (sec)") 
ylabel ("SOC(%)") 
xlim([30000 30100]) 


yticks ([ -0.046 -0.044 -0.042 -0.04 -0.038 - 
0.036 -0.034]) 

yticklabels({ '75.4', '75.6', '75.8', '76', '76.2', 
76.4' , '76.6' }) % set the label to percentage 


%%%Generator output graphs 

figure (3) 
plot(Igen) 

title( "Generator Current ") 
xlabel("Time (sec)") 
ylabel ("Current (A)") 
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%%%Super Capacitor plotting 

figure (4) 
subplot (2,1,1) 
plot (ISC) 

title ("Super Capacitor Current(A) ") 

xlabel("Time (sec)") 
ylabel ("Current (A)") 

subplot (2,1,2) 

plot(SOC_UC) 


title ("Super Capacitor SOC ") 
xlabel("Time (sec)") 
ylabel ("SOC(%)") 

figure (55) 

SOC_UC_Inv = S0C_UC*-1; % inverse the plot 
plot(SOC_UC_Inv) 

o 

o 

yticks([-0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 
- 0.1 0 ]) 

yticklabels ({'O', '20', '30', '40', '50', '60', '70', 

'80' , '90' , '100' }) % set the label to percentage 

title ("Super Capacitor SOC ") 
xlabel("Time (sec)") 
ylabel ("SOC(%)") 

% % invert the graph , enahnce view 

SOC_UC_Inv = SOC_UC*-l; % inverse the plot 

plot(SOC_UC_Inv) 

title ("Super Capacitor SOC ") 

xlabel("Time (sec)") 

ylabel ("SOC(%)") 

xlim([30000 30100]) 
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yticks([-0.35 -0.3 -0.25 -0.2 -0.15 -0.1]) 
yticklabels ({'65.5', '75', '75.5', '80', '85.5', '90 

' }) % set the label to percentage 

title ("Super Capacitor SOC ") 
xlabel("Time (sec)") 
ylabel ("SOC(%)") 


%% Comparison 

%%ISC and IBatt 

figure (5) 
plot (ISC) 
hold on 
plot(Ibatt) 
hold off 

legend( "Isc" , "Ibattt" ) 
title ("Ibatt and ISC") 
xlabel ("Time (sec) ") 
ylabel ("Current (A) ") 

%%%% Current comparisons 

figure (7) 

plot(w(1,:), "r" ) % wave curent 

hold on 

plot(Ibatt, "b" ) %battery current 
plot(Igen ,"g") %Generator current 
plot( ISC,"c") % SC current 

plot(Igen+(Ibatt) + (ISC), "k" ) % total current 

legend( "I wave" , "Ibatt ", "Igen" , "ISC" , "Total 
Current" ) 

title (" Currents Combine plot") 
hold off 

xlabel ("Time (sec) ") 
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ylabel ("Current (A) ") 

figure (30) 
plot(w(1, :), "r" ) 
hold on 

plot(Ibatt, "b" ) 
plot(Igen , "g" ) 
plot( ISC, "c" ) 

plot(Igen+(Ibatt)+(ISC), "k") 
xlim([3400 3500]) 

legend( "I_{wave}" , "I_{batt}" , "I_{gen}" , "I_{SC}" 

, "Total Current") 

title (" Currents Combine plot") 

xlabel ("Time (sec) ") 
ylabel ("Current (A) ") 
hold off 

figure (31) 
subplot (3,1, 1) 
plot(Igen) 

title( "Generator Current (A) ") 

xlabel( "Time(Sec)" ) 
ylabel ("Current (A) ") 

%ylim([-250 250]) 
xlim([3400 3500]) 


subplot(3,1,2) 
plot(ISC) 

title ("Super Capacitor Current(A) ") 

xlabel ("Time (Sec) ") 
ylabel ("Current (A) ") 

%ylim([-250 250] ) 
xlim([3400 3500]) 

subplot(3,1,3) 
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plot(Ibatt) 

title( "Battery Current(A) ") 

xlabel ("Time (Sec) ") 
ylabel ("Current (A) ") 

%ylim([-250 250]) 

xlim([3400 3500]) 

Geny_max = max((Igen(3400:3500) ) ) 
Batt_max = max((Ibatt(3400:3500))) 
SC_max = max((ISC(3400:3500) ) ) 
hold off 

B. SHIPPOWER FUNCTION 

function [Pe,v,DSP] = shippower(v_kn) 

% % MF Oppedal 
% % 

% % numbers of the battery, add SC and 
simulink, 

% Drag calculation for MF Oppedal 


%Standard Parameters 

rho_h2o = 1025; %kg/m A 3 density of seawater 
rho_air = 1.22; %kg/m A 3 desity of air 
g = 9.81; %ms A -2 

vis = 1.15*10 A -6; %ITTC Fresh Water and 
Seawater Properties avg viscosity (kin. 
vis./rho) 

% v_kn = 15; %1inspace(1,15,30); 
v = v_kn * 0.514444; %conv. to ms A -l 

%Dimensions in m 
B = 16 ; %breadth 
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Lpp = 105 ; %m perpendicular length source: 
Kullmann 

LWL = Lpp * 0.97 ; % 97% of Lpp, according to 
literature 

DSP = 1705*10 A 3 ; % kg Displacement (weight) 
PLD = 600*10 A 3 ; %kg payload max. 
nabla = (DSP+PLD)/rho_h2o ; % Displacemed 
volume on maximum cap. 

A_m = 18*2.3 ; % assumed wetted cross section 

area 

T = 3.3 ; % m from www.marinetraffic.com 

(average value) 

A_x = 18*10 ; % aussumed cross section area 
above water 

Re = (5e8);%(v * LWL)/vis % Reynold No. 

Fn = v/( (g*Lpp) A 0.5) ; %Froude number 

S = 0.99 *( (nabla/T)+1.9*LWL*T) ;% wetted 
surface area Mumford 


% %%%%%%RESISTANCE COEFFICIENTS%%%%%% 

C_B = nabla /(LWL*B*T); % block coeff 

C_F = 0.075 / ((logl0(Re) - 2 ) A 2); 

;%frictional Resistance 

C_X =0.8; % for bulky geometries —> acc. to 
ITTC 

C_AA= C_X * ((rho_air*A_x)/(rho_h2o*S)); % air 

Res . 

C_P = nabla/(LWL*A_m); %prismatic coefficient 
% residual coeff. 

C_R = 1* 10 A -4; %from Harvald/Guldhammer Chart 
% Form Factor k 

k = -0.095 + 25.6 *(C_B/(LWL/B) A 2 * (B/T) A 0.5) 

%Length Displ. Ratio 

LDR = Lpp/(nabla A (1/3) ); %molto importante to 
determine which diagram to use 
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C_RA = 0.044 *((((150*10 A -6)/LWL) A (1 / 3 ) ) - 

10*(Re A (-1/3)))+0.000125; % roughness allowance 

%TOTAL COEFF 

C_T = (1+k)*C_F + C_R + C_AA + C_RA; %Total 

Res . 

R_T = 0.5 * C_T * rho_h2o * S * v A 2; 

Pe = (R_T* v)/0.51; %0.51 efficiency from mech 

to electrical 

end 


C. POWERSCRIPT FUNCTION 

% %%% MODEL 4_based on current for Battery and 
SuperCap 

% % Regulation with sinusoidal exosystem. 

% % x=[Q_tilde, Igen, Ibatt, SOCbatt-SOCrated, 
Isc,SOC_SC]' 

% % u=[DeltaPgen, Deltalbatt, Deltalsc]' 

% %%% STATE SPACE MODEL 
function 

[e, x, u, w] =P owerscript(I_wave,wl,w2,theta0,FS,Ig 
en,t,Voc,Q_batt,Vmax,SC_Cap,T_wave); 

Ts = 1/FS 

gamma = -(T_wave/(3600*200)) % Large gamma = 

battery size smaller %200 

beta = -(T_wave / (Vmax*45*7)) % , S C = 

0.00001 , 5 parallel branches 

F = [0.8, 1, 1,0,1,0 ; 

0 , 1 , 0 , 0 , 0 , 0 ; 

0, 0, 1, 0, 0, 0 ; 

0, 0, gamma, 1,0,0; 

0, 0, 0, 0, 1, 0 
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0 , 0, 0, 0,beta,1]; 

G = [ 0,0,0; 

1 , 0 , 0 ; 

0 , 1 , 0 ; 

0 , 0 , 0 ; 

0 , 0 , 1 ; 

0,0,0] ; 

B = [-1, 0; 

0 , 0 ; 

0 , 0 ; 

0 , 0 ; 

0 , 0 ; 

0 , 0 ]; 

H = [0, 1, 1,0,1,0]; 

J = [ 0,0,0] ; 

D = [-1,0]; 

% EXOSYSTEM 

A= [cos(thetaO) , -sin(thetaO); sin(thetaO) 
cos (thetaO)]; 


ooooooooooooooooooooooooooooooooooooooooo 


% dimensions 

Nx=size(F,1); 
Nu=size(G, 2) ; 
Ny=size(H,1); 
Nw=size(A,1) ; 


% state size 
% input size 
% output size 
% exosystem state size 


%%% Cost Parameters 

%the higher the cost means, more peantly, 
use 


D, O 
O O 


less 
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rg=l; rbatt=l; rsc=l;% less weight on the rate 
of change for the delta values as we want them 
to go to steady state quicker 

qxv=l; % weight on the error (velocity) higher 
= lower ripple in the error 
qgen = 50; % weight on generator 

qlbatt =50;% weight on battery %the battery has 
give out more power as compared to geny % 
inverse to qsoc 

qSOC_batt=l; % weight on Battery SOC % 

qISC = 1; % 
qSOC_SC =1; 

% Since |e| A 2 = z'C'Cz+v'D'Dv+2z'C ' Dv then 

diag_q = [qxv qgen qlbatt qSOC_batt qISC 
qSOC_SC]; 

% Q=qe*H' *H; 

Q=diag(diag_q)+(H'*H); 

R=diag([rg,rbatt,rsc])+qxv*J' * J; 

E=qxv*H'* J; 

% Riccati Equation 

[K,P,CLP]=dlqr(F, G, Q, R, E) ; 


%x (t) 

Igen_max = inf%100 %100;% bounds on Igen no 
bounds on the generators , so to get a quicker 
system reponses as well as to see the effect on 
the rae of change of the gen 
Igen_min = inf%100; %bounds on Igen 

%Battery constraints 

% assuming the Battery starts at 50%, 50% = 0 
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%rate of change for SOC for the batteries.. 
time to reach steady state , 

DSOC_batt_max= 0.8 ; %max 90% 

DSOC_batt_min= 0;%-0.2;% min 20% 

Ibatt_max = 800 %400 % mac discharge is 400A 
Ibatt_min = 800 

%Super_cap Constraints 

ISC_max = 240*4; 

ISC_min = 240*4; 

% assuming the Supercapacitor starts at 50%, 
50% = 0 

SC_SOC_max = 0.9; % 

SC_SOC_min = -0.1;% 

% constraints on output XV 

Xv_min = 20; % constraints on the error in 
0 tilde, speed vary 
XV_max = 20; 

Ax=[0,0,1,0,0,0; % IBattery 
0,0,-1,0,0,0; % IBattery 

0,1,0,0,0,0; % Gen 
0, -1,0,0,0,0; % Gen 
0,0,0,1,0,0; % Battery SOC 
0,0,0,-1,0,0; % Battery SOC 
0,0,0,0,1,0; % ISuperCap 
0,0,0,0,-1,0; % ISuperCap 
0,0,0,0,0,1; % SuperCap_SOC 
0,0,0,0,0,-1; % SuperCap_SOC 
1,0,0,0,0,0; % XV 
-1,0,0,0,0,0]; % XV 


70 


bx=[Ibatt_max;Ibatt_min;Igen_max;Igen_min;DSOC_ 
batt_max;DSOC_batt_min;ISC_max;ISC_min;SC_SOC_m 
ax;S C_SOC_min;XV_max;Xv_min]; % 

% Constraints on Variables 
% u (t) 

DPgen_max=50 ; % max on |DeltaPgen| % should be 
the rate of change for delta Pgen .. Smaller 
Dpgen Max, slower rate of change for the 
generator 

% the lower the value for DPgen, the slower it 
will reach a steady state value 

DPbatt= 500; 

DPSC = 5000; 

Au=[1,0,0; % Constraint on the Generator rate 
on change of 
- 1 , 0 , 0 ; 

0 , 1 , 0 ; 

0 ,- 1 , 0 ; 

0 , 0 , 1 ; 

0,0,—1]; % Constraint on the Generator rate on 
change of 

bu=[DPgen_max;DPgen_max;DPbatt;DPbatt;DPSC;DPSC 

]; 

N=20; % length of Prediction Horizon, 2 sample 
wave, based on the samping Freqency 

[Linv,Fmpc,Gu,Hu,hu, options, 

iAO]=mpcsetup(N,F,G,Q,R,E,P,Ax,bx,Au,bu); 

% Solve Francis Equation 
% T is Nx by Nw, L is Nu by Nw 

[T,L]=Francis(F,G,H,J,A,B,D,Nx,Nu,Ny,Nw); 
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% Check Francis Equation: these terms should be 
zero 

F*T+G*L+B-T*A; 

H*T+J*L+D; 

% Constraints on v=u-Lw are of the form 
Gu*V<Hu*z(t)+hu+Hw*w (t) 

% Just follow the notes: 

Hw=[]; 
for k=l:N 
Hw=[Hw;-Ax*T*A A k] ; 
end 

for k=l:N 
Hw=[Hw;-Au*L*A A k] ; 
end 

S-2-2-2-3-2-2-3-3-2-2-2-2-2-2-2-2-2-2-2-3-2-2-2-2-2-2-2-2- 

ooooooooooooooooooooooooooooo 

%%% MPC Control Simulation 

Nsim=length(t)+1; % length of simulation 

n=l:Nsim; 

% Initial State 

x(:,1)=[0;wl(1);0;0;0;0]; % Q tilde, Pgen, 
IBatt, SOC_B, ISC,SOC_SC 
w ( :,1) = [wl(1) ;w2(1) ] ; 
z ( :, 1) =x ( : , 1) -T*w ( :, 1) ; 

for t=l:Nsim-l; 

[V, status, iA, lambda] = mpcqpsolver(Linv, 
Fmpc'*z(:,t), -Gu, -Hu*z(:,t)-hu-Hw*w(:,t), 
double.empty(0,1) , double.empty(0,1), 
iAO,options); 
v(:,t)=V(l:Nu); % MPC 

e (t)=H* z(:, t)+J*v(:,t); 
z(: , t+ 1)=F* z(:,t)+G*v(:,t); 
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w(:,t+1)=A*w( :, t); 

% back to riginal x and u 
x(: , t)=z(:,t)+T*w (:,t); 
u(:,t)=v(:,t)+L*w ( : , t) ; 
end 

end % End of Function 

D. MPCSETUP FUNCTION 

function [Linv,F,Gu,Hu,hu, options, 

iAO]=mpcsetup(N,A,B,Q,R,E,P, Ax,bx,Au,bu) 

% [Linv,F,Gu,Hu,hu, options, 

iAO ] =mpcsetup (N, A, B, Q, R, E, Ax, bx, Au, bu) 

% by R. Cristi, 8/1/2019 
% INPUTS: 

% N=length of time horizon 
% A,B = state equations 

% Q,R,E= quadratic penalties for state,input 
and crossterm (2x'Eu) 

% P = last state penalty from ARE 
% A_x, b_x, A_u, b_u constraints on x and u 
% OUTPUTS: 

% Optimization of U'HU+2x(t)'FU, w.r.t. U 
% Linv=a Choleski representation of H 
% Constraints of the form 
% Gu*U<Hu*x(t)+hu 

% options and iAO parameters for the algorithm, 
given as default. 

% At each iteration "t" call the algorithm 
% [U, status, iA, lambda] = mpcqpsolver(Linv, 

F'*x(:,t), -Gu, -Hu*x(:,t)-hu, 
double.empty(0,1) , double.empty(0,1), 
iAO,options); 

% u(:,t)=U(1:Nu); % Nu=size(B,2), number of 
inputs 

e-s-s-e-e-s-s-e-e-s-s-e-e-s-s-e-e-s-s-e-e-s-s-e-e-s-s-e-e-s-s-e-e-s-s-e-e-s- 

00000000000000000000000000000000000000 
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dimension of the state 
number of inputs 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

oooooooooooooooo 


Nx=size(A,1); % 

Nu=size(B,2) ; % 


%%%% Matrix Sx 

Sx=eye(Nx); 
for k=l:N 
Sx=[Sx; A A k]; 
end 

%%%%%% Matrix Su 

Suk=zeros(Nx, Nu*N); % first block 

Su=Suk; 

for k=l:N 

Suk=[(A A (k-l))*B, Suk ( :,1 : (end-Nu)) ]; 
Su=[Su;Suk]; 
end 


% Matrices Qbar, Rbar, Ebar 

% Note: use "kronecker" product to form block 
diagonal matrices: 

Qbar=kron(diag([ones(1,N),0]), 

Q)+kron(diag([zeros(1,N), 1]), P); 

Rbar=kron(eye(N) , R) ; 

Ebar=kron(eye(N) , E) ; 

% Sxbar, Subar in the notes: 

SufirstN=Su(1:N*Nx, :); SxfirstN=Sx(1:N*Nx,:); 

% first N locks 

% Sx, Su underlined in the notes: 

SulastN=Su(Nx+1:end, :); 

SxlastN=Sx(Nx+1:end,:); % last N blocks 
% Matrices H,F,Y 
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H=Su'*Qbar*Su+Rbar+SufirstN'*Ebar+Ebar'*Sufirst 

N; 

F=Sx'*Qbar*Su+SxfirstN'*Ebar; 

Y=Sx'*Qbar*Sx; 

% Setup Constraints of the form 
% Gx*z<Hx*x+hx 
% Gu*z<Hu*x+hu 

% SxlastN, SulastN are Sx and Su without the 
first block (ie the first Nx rown) 

Axbar=kron(eye(N) , Ax) ; 
bxbar=kron(ones(N,1) , bx) ; 

Aubar=kron(eye(N), Au); 
bubar=kron(ones(N, 1) , bu) ; 

% constraints of input U from state x: 

% Gx*U<=Hx*X+hx 

Gx=Axbar^SulastN; 

Hx=-Axbar* SxlastN; 
hx=bxbar; 

% total constraints on U: 

% append constraints Aubar*U<bubar to 
constraints from X: 

% Gu*U<Hu*X+hu; 

Gu=[Gx;Aubar]; 

Hu=[Hx; zeros(size(Aubar,1), Nx)]; 
hu=[hx;bubar]; 

%% Parameters for mpcQPsolver 
Ng=size(Gu, 1) ; % number of constraints 

% Choleski decomposition for H 
L = chol(H, 'lower'); 

Linv = L\eye(size (H,1)) ; 
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% SET options for "mpqpcsolver" (in 2019b, just 

copy from default in help file): 

options.DataType= 'double' ; 

options.Maxlter=2 0 0; 

options.IntegrityChecks=true; 

options.FeasibilityTol=1.0e-06; 

options.Integrity.Checks=true; 

iA0=false(Ng, 1) ; 

% END 

E. FRANCIS FUNCTION 

function 

[T, L]=Francis(F,G,H,J,A,B,D,nx,nu,ny,nw) 

% This function returns the solution of 
% the discrete time Francis equations of 
% linear regulation. 

% F*T+G*L+B=T*A 
% H*T+J*L+D=0 
% 

% The plant state is X which is NX dimensional. 
% The plant control is U which is NU 
dimensional. 

% The plant output is Y which is NY 
dimensional. 

% The exosystem state is X_ which is NW 
dimensional. 

% The plant dynamics is 
% $$ x A +=F*x+G*u+B*w $$ 

% The plant output is 
% $$ y=H*x+J*u+D*w $$ 

% The exosystem dynamics is 
% $$ w A +=A*w $$ 

% If there is no unique solution to the Francis 
equations 

% then a least squares solution is found. 

% The result will be unsatisfactory 
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% 


if NY>NU. 


% 


% We use the identity 
% (A*Z*B)rs= Zrs*kron(A',B) 

% where A, Z, B are matrices of dimensions 
% compatible for multiplication 
% and rs means row stacked. 

templ=kron([F, G; H, J] . ',eye (nw) ); 
temp2=templ- 

kron ( [eye(nx),zeros(nx,nu) ;zeros (ny,nx+nu) ] 

) ; 

temp=-reshape([B; D] . ' , 1 , (nx+ny)*nw); 
temp=temp*pinv(temp2); 

T=reshape(temp(1,1:nx*nw) , nw, nx) . ' ; 

L=reshape(temp(1,nx*nw+l: (nx+nu)*nw) ,nw,nu) 


% 


Copyright A. J. Krener 2019, 
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